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We discuss algorithms applicable to the numerical solution of 
second-order ordinary differential equations by finite-differences. We 
make particular reference to the solution of the dissipative particle dy- 
namics fluid model, and present extensive results comparing one of 
the algorithms discussed with the standard method of solution. These 
results show the successful modeling of phase separation and surface 
tension in a binary immiscible fluid mixture. 

I. INTRODUCTION 

Computer simulations are a test-bed for theory, permitting experiments over 
which we have nearly complete control. The effectiveness of these computer ex- 
periments is limited by the power of available computers, and by their inherently 
discrete nature since the physical and chemical phenomena we study are continu- 
ous. Fortunately, mathematics has provided us with a framework into which we 
can place this problem: finite-differences and their continuous analogy, differential 
equations. 

In the present paper, we shall be concerned with algorithms for solving second- 
order ordinary differential equations of the form 

x = f(x, x), (1) 

where x denotes a scalar or vector position, and dots indicate derivatives with re- 
spect to time. In particular, we are concerned with determining the motion of a 
system of particles moving. according to the dissipative particle dynamics (DPD) 
computational fluid modeltl (see Section |l|). The large number of equations and 
their non-linear nature make direct analytic mathematical solution infeasiblc, leav- 
ing numerical solution by finite-differences as the most appropriate approach. The 
forces in our computational fluid model are unusual in that they have a stochas- 
tic component, which depends on the size of the timestep in our finite-difference 
algorithm. 

Traditional methods for solving Eq. (Q) can be found in textbooks on numerical 
analysis ,mj and in .papers and texts on molecular dynamics and related compu- 
tational methods. ETEj Algorithms from numerical analysis textbooks are usually 
unsuitable for our purpose because of the unstable nature of the motion of a many- 
body system, or an assumption that the forces do not depend on the timestep. Most 
of the other algorithms we find assume conservative forces, i.e. x — f (x). 

Throughout this paper, we use the notation x n = x(t) and x n +i — x(t + h) for 
the value of a variable (vector or scalar) at successive timesteps, where h > Q is the 



timestep width. Variables with tildes denote temporary quantities, vectors are in 
bold, and matrices are in a sans serif font. 



II. DISSIPATIVE PARTICLE DYNAMICS 

Hoogerbrugge and Koelman proposed DPD0 as a novel particulate model for the 
simulation of complex fluid behavior. This model is particularly well-suited to model 
multiphase flows, flow-in porous media, colloidal suspensions JHia microemulsions, 
and polymeric fiuids.tU The traditional approach of continuum fluid mechanics has 
met with limited success, and so many new micro- and mesoscopic approaches 
have been considered. In principle, molecular dynamics (MD) is the most accurate 
microscopic approach, although in practice it is too slow in both its quantum (Car- 
Parinello) and classical forms because of its excessive detail. Discrete methods 
developed from lattice-gas automata (LGA) have had-some success, but they too 
have problems, such as lacking Galilean invariance.EilO 

DPD was developed in an attempt to capture the best aspects of MD and LGA. 
It avoids the lattice-based problems of LGA, yet maintains an elegant simplicity and 
larger scale that keeps the model much faster than MD. This simplicity also makes 
DPD highly extensible, such as for including the interactions of complex molecules 
or modeling flow in an arbitrary number of spatial dimensions. The key features of 
the basic model are that the fluid is grouped into packets, termed "particles" , and 
that mass and momentum are conserved but energy is not. Particle positions and 
momenta are real variables, and are not restricted to a grid. 

There are two sequential steps to the action of the original model proposed by 
Hoogerbrugge and Koelman:El (i) an infinitcsimally-short impulse step 

Pz,n+1 — Pz,n H~ ^ ^ ^ij&ijj (2) 

and (ii) a propagation step taking time h 

x i.n+l = x i.n Pi.n+1; (3) 

m 

where Pi and are the momentum and position of particle i, is the unit vector 
pointing from particle j to particle i, m is the mass of each particle, and 

3m (l - ^fj 

n tj = - — — [Uij -ufa- pj) ■ eij] (4) 

Til c p 

within the cut-off radius r c . The quantity denotes the distance separating par- 
ticles i and j, and p is the mass density; note also that the normalization given is 
correct only in two dimensions. The random variable TLij represents the conserva- 
tive and stochastic effect of the collisions and gives rise to fluid pressure, while the 
second, dissipative, term inside the square brackets yields fluid viscosity. 

Espanol and Warren's analysis!!] showed that the original DPD model does not 
satisfy detailed balance, so the equilibrium states (if they exist) cannot be simply 
characterized. Detailed balance is the condition equating the rates of forward and 
backward transition probabilities in a dynamical system, and is a sufficient (but 
not necessary) condition guaranteeing that the system has a (Gibbsian) equilibrium 
stateO'ta Espanol and Warren formulated a Fokker-Planck equation and equivalent 
set of stochastic differential equations which lead to a similar model, 
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< i¥* (5) 

In these equations, p^, x,, and m, denote the momentum, position, and mass of 
particle i, and Ffj is a conservative force acting between particles i and j while 
and are the dissipative and random forces; dWij = dWji are independent 
increments of a Wiener process. By Ito calculus 

dW tj dW k i = (SikSji + 5u5jk) dt, (6) 

and so dWij is an infinitesimal of order | and we can write dW^_= Oij^fdi, where 
9ij = 9ji is a random variable with zero mean and unit varianceO With an appro- 
priate choice for the form of therferces we find detailed balance is satisfied by this 
continuous-time version of DPD,li3 and so equilibrium states are guaranteed to ex- 
ist and be Gibbsian. To ensure that the associated fluctuation-dissipation theorem 
holds, the forces assume the following forms 

Fg = auijeij, (7) 



F S = -^4(^- v y)^> ( g ) 

and 

Fg = owfc-ey, (9) 

where = Pi/rrij — Pj/rrij is the difference in velocities of particles j and i, is 
the unit vector pointing from particle j to particle i, and u)ij is a weighting function 
depending only on the distance separating particles i and j. The constants a, 7, 
and a are chosen to reflect the relative importance of the conservative, dissipative 
(viscous), and random components in the fluid of interest. As a consequence of 
detailed balance and the fluctuation-dissipation theorem, 7 and a are related to 
Boltzmann's constant kg and the equilibrium temperature T by 

2 

— = 2k B T. (10) 

7 

In order to remain as close as possible to the original DPD model, we choose the 
friction weight function to be 

w« = i-Ji (ii) 

<c 

within the constant cutoff length r c > 0, where is the distance between particles 
i and j. Adding Eqs. (@)-(|), 

the total force is 



a-JUij {eij -Vij) H 



(12) 



In summary, the main changes from the original DPD model are the specification 
of the motion in terms of differential equations instead of an impulse step followed 
by coasting; the separation of the conservative, dissipative, and random forces, 
with the strength of each being controlled by the new constants a, 7, and a; the 
insertion of an extra factor of w,j in front of the dissipative force; the specification of 
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the thermodynamic temperature in terms of the model parameters; and the square 
root dependence of the random force on size of timestep. 

In order to model binary immiscible fluids, we adopt the simplest approach of 
introducing a new variable, the "color", by analogy with Rothman-KellerB3 When 
two particles of different color interact we increase the conservative force, thereby 
increasing the repulsion. That is, 

{«o if particles i and j are the same color 
a\ if particles i and j are different colors, ^ ' 

where «o and a± are constants with < cto < cki- As for the single-phase DPD 
fluid, the Navier-Stokes equations are obeyed within regions of homogeneity in each 
of the two immiscible, fluids, while detailed balance is preserved, at least in the limit 
of continuous time.Eil We would like to emphasize that with the sole exception of 
the choice of finite-difference algorithm, this model is identical to that used in our 
most-recently published simulationsHj 



III. EVALUATION CRITERIA 



A good finite-difference method should be fast, stable, accurate, easy to imple- 
ment, and require little storage. Since for most large systems the force evaluation 
dominates the computation time, the fastest algorithm is the one that evaluates 
the force the fewest times for a given accuracy. In general, this means that stable 
methods which allow large time steps are preferable. In this context we must be 
careful to distinguish the stability of a method from the stability of the system 
being studied, as the dynamical systems we wisk-jta, study are extremely unstable, 
with small perturbations growing exponentiallyJ§BE3 A stable method will respond 
better to large sizes of timestep than an unstable method. The accuracy of an al- 
gorithm is difficult to measure in practice because of the extreme instability of our 
dynamical systems. Practical accuracy is best quantified by requiring a method to 
give a trajectory in phase space that is physical for the system (i.e. similar to that 
observed in experiments), and for which calculated properties of the system are 
close to their theoretical values. An idea of the accuracy of a method can also be 
obtained by considering the local truncation error, expressed in terms of the order 
of the timestep occurring in the first truncated term of the corresponding Taylor 
series. Because the forces considered are not smooth, low order algorithms may 
possibly be more accurate than high order algorithms. 

What jCpiteria can we use to evaluate the suitability of a method? The 
literature3i§E2l claims that a good finite-difference method will be area preserving, 
time reversible, energy conserving, and ergodic. The first three criteria are closely 
related. We need to be careful in distinguishing between the properties of the finite- 
difference algorithm and those of the system it is solving. For example, when we 
describe an algorithm as being time reversible, we are not saying that it corrects 
or compensates for the lack of time symmetry in a non-conservative system, but 
rather that when applied to a conservative system this finite-difference algorithm 
is able to retrace the precise sequence of steps in the backward evolution which it 
traversed originally in the forward time direction. A similar confusion may arise 
when we consider ergodicity. Is this not a property normally associated with the 
rules of interaction of systems? What does it mean to say that a finite-difference 
algorithm is ergodic? When we describe an algorithm as being ergodic, we mean 
that an ergodic system updated by this algorithm will still sample the full phase 
space; a non-ergodic algorithm updating the same system will not. 
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We also need to consider the relevance of these four criteria to non-conservative 
forces. For example, is there any use in having an algorithm which conserves energy 
when our DPD interactions conserve only momentum? We argue that unnecessary 
losses are indeed undesirable; a finite-difference algorithm which loses energy by con- 
sistently underestimating velocity would exhibit lower temperature than a similar 
algorithm which does not. 

An algorithm which preserves the volume of phase space will have a unit Jacobian 
for the transformation of variables from one timestep to the next. Time reversibility 
can also be checked analytically (ignoring for the moment the properties of the 
forces). Energy conservation normally follows from time reversibility, but it can 
also be verified for a conservative system, such as the (one-dimensional) simple 
harmonic oscillator 

x = -x, (14) 

where the total energy is E = (x 2 + x 2 )/2. 

Ergodicity can be checked by calculating the equilibrium temperature and pres- 
sure of an homogeneous fluid, and comparing the velocity distribution with a normal 
(Gaussian) distributions As this is not a sufficient condition, in the present paper 
we also simulate a binary immiscible fluid using each of the algorithms, testing 
Laplace's law and comparing the rate of growth of domains with the extensive ex- 
perimental and theoretical literature, and with our earlier simulations! 22 '! 24 ! using 



the standard DPD finite-difference algorithm (see Section IV A). It has also been 
noted in DPD simulations of an ideal gas (i.e. a = 0) using this algorithm that the 
radial pair-correlation function differs from the theoretical result Hj A good finite- 
difference algorithm would remove this inconsistency, so we verify this as well. 



IV. ALGORITHMS CONSIDERED FOR CONTINUOUS-TIME DPD 

For molecular dynamics simulations, the literatur|er*eepmmends either one of the 
Verlet schemes or a multi- value predictor-corrector JaflE2l Runge-Kutta and extrap- 
olation methods are considered unsuitable because they typically require a large 
number of force evaluations per timestep, which is inefficient for a many-particle 
system. For the same reason, methods which consider the second order differential 
equation directly are usually preferred to methods which break it into two first order 
equations. 

We do not consider implicit methods (in which a variable appears on both sides 
of a single equation) since they are computationally infeasible for any system large 
enough to be of interest. This is so because the simplest implicit method for a 
system of N particles -requires solution of N interdependent linear equations tak- 
ing 0(iV 3 ) operations^! compared with the 0(N) operations of an explicit method. 
The need to solve these interdependent equations also makes implicit methods much 
more complicated than explicit methods. However, where the differential equations 
are very stiff (have several greatly differing time scales) the greatly-increased algo- 
rithmic stability could outweigh the disadvantages. Having had reasonable results 
with explicit methods, we do not believe our system is so stiff as to require implicit 
schemes. 

We have obtained good results by starting our simulations with all particles at 
rest, with the initial position of each particle typically selected from a random 
variable uniformly distributed across the simulation space. Particles are assumed 
to have been at rest for some time, so prior positions are the same as initial positions, 
and prior velocities and higher derivatives are all zero. In our notation x n = x{t) 
and x n +i — x(t + h), the initial conditions at t = correspond to n = while 
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negative values of n represent values prior to the start of the simulation. To give an 
idea of the relative accuracy of each algorithm we report the local truncation error 
in position as the order of the timestep appearing in the first term omitted from 
the Taylor series. 



A. Basic methods 

The original algorithm proposed by Hoogerbrugge and KoclmarE is a modified 
Euler scheme 

.1 n 

Xn+l 

As much researpiui j— .J— ||Jjas-. been published using 

this area-preserving scheme Jd'0£§E3EirE3 it is important to compare it with the 
other methods. The local truncation error is 0(h 2 ). 



/ (*£n; %n j ft) 
*^n "t - hx n 

— *^n hXyi-\-\ 

— x n + h <x n + hx n ) . (15) 



B. Verlet-based methods 

Due to the velocity-dependent nature of the force (Eq. (|l^)), the traditional 0(/i 3 ) 
time-reversible velocity- Verlet algorithm 

hi- h " 

%n-\-l %n "T ht I X n T" ^ 

= / (#n+l) 

X n +l = %n + ^ + *n+l) (16) 

cannot be used. Ferrarioi has suggested two modified schemes to get around this 
problem. The first is 

Xn f (*Cfis %n i tl) 
<^n— 1 

Xn ~ 2h 

X n +1 = X n -1 + 2hx n (17) 



and the second 



hi- h- 

±n+l = in + ^ (in + X n +\) ■ (18) 



These two Q(/t 3 j methods are hereafter referred tOr-as Ferrario I (Eq. (|i7|)) and 
Ferrario II (Eq. (|18[)) respectively. Groot and Warrencl recently proposed a different 
velocity estimate for use with the velocity- Verlet algorithm. Their method is 



; . , . h.. 
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X 

Xn+1 
Xn+l 

This is conceptually similar to the leap-frog form of the Verlet algorithm, for if 
A ^ 1 the force is evaluated when the positions and velocities are out of phase. A 
reasonable range for A would be from to 1; we consider only A = | and A = 1. 
For conservative forces Eq. ( |l9| ) reduces to Eq. (|l6|); otherwise with A = 1 its local 
truncation error is 0(h 3 ), and 0(h 2 ) for A ^ 1. 



- — - <3j/), I hhx fi 

= f (x n +l, X, h) 

= x n + - (x n + x n+ i) . (19) 



C. Runge-Kutta methods 



In order to confirm the poor performance of Runge-Kutta schemes, we,-consider 
a two-stage explicit Runge-Kutta method of the form proposed by Geaia for the 
direct solution of second order differential equations: 

fi / ^£tu x n , —h 

2 u( ■ h i 

x = x„ + -h \x n + - fi 
2 ~ 

X = X n + -kfl 

1. 



f 2 = f [ x, x, -h 



Xn+l =X n + h ^X n + J (/l + h^j j 

(A + 3/ 2 ) . (20) 



h 

Xn+l — Xn T" 



In applying this 0(/i 4 ) method to our system, we must be careful to choose the 
appropriate length of timestep for each force evaluation. 



D. Multi- value predictor-correctors 



Gear has proposed, a notation for the general class of multi-value predictor- 
corrector algorithms The variables to be stored from step to step are kept in 
a column vector, usually in either the Nordsieck (N) representation 

yn(A0=(zn, hi n , yl n , ^X^ > , ^M, . . .) (21) 

or the force (F) representation 

Yn \-T ) — I X n j i^X n , ^ Xji , ^ Xn— 1 j ^ Xn — 2 )••*!• V / 

The iV-representation is more convenient for changing timestep and order, while the 
^-representation simplifies calculation. Transformation from one representation to 
the other is straightforward, involving a multiplication by a transformation matrix 
T. Gear claims that transformation from one representation to another does not 
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affect the truncation error or stability of the method, but does affect the round-off 
error, and amount of computation and storageE'EH 

To take a step forward in time using a predictor-corrector method, we first mul- 
tiply the column vector y n -i with a matrix A in an explicit prediction step 

y n ,o=Ay„_i > (23) 

and then the forces are evaluated at the positions and velocities in this vector, 

A corrected vector is obtained by adding a multiple of the difference between the 
predicted and the calculated force to the existing vector 



m '■ yn,(m-fl) 







(25) 



It is possible to iterate the force evaluation (E) and correction (C) steps after the 
initial prediction (P), leading to a P(EC) q or P(EC) q E algorithm. The subscript 
to denotes each evaluation-correction step from the first (to = 0) through the last 
(to = q — 1, i.e. < to < q). It can be shown that in the. limit of a large number of 
iterations (g> 1), these algorithms are time-reversible 

The predictor matrix A is usually chosen to be the Pascal triangle matrix 



A(N) 



Z 1 


1 


1 


1 





1 


2 


3 








1 


3 











1 



V 



(26) 



(in the N- representation), which means that Eq. ( f23[ ) predicts according to the 
familiar Taylor series expansion. The column vector 1 is chosen by ac&ttracy and 
stability arguments. The corrector vectors for 3- and 4- value methodsEa to solve 
second order differential equations in the Af-representation are 



W) = ( 3, 1, 1 



and 



(27) 



1(A0 



1 5 1 

6' 6' ' 3 



(28) 



giving local truncation error of 0(/i 4 ) and 0(h 5 ) respectively. The correction vectors 
are identical in the ^-representation except that all components beyond the third 
are zero. BO 

We would like to consider these two methods further modified by the insertion of 
a factor A (typically ranging from to 1) into the velocity prediction: 



A(AT) = 



1111 
1 2A 3 
13 
1 



(29) 



(this matrix being for the 4- value method). While this predicted value of velocity 
will be used for the first force evaluation, the Taylor-predicted velocity will be used 
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in the correction step and subsequent iterations so that the methods will benefit 
from the modified velocity estimate and yet will still converge to the correct solution. 
It can be shown that both Eqs. (p6| ) and ( p9| ) are area-preserving transformations. 
Written explicitly in the P(EC) ''F-representation, we have 



P : 



x n -i + h ( i„_i + 2*™-i 



E m • ^n,(m+l) 



/ (x n;0 , X71-1 + \hx n -i, h) if m = 



if m > 



and 



C ,i 



^.(m+l) — %n,m H~ g (^n,(m+l) 



(30) 



./■„.„ = x n -x + h [ x n -i + - (4aS n _i - £„_ 2 ) 



■ r ' : <I i„ i0 = i n _i + - (3x„_i - x n _ 2 ) 

{ %(n-l),0 — x n-l 



f ( a^n.o, in-i + A- (3i„_i - £„_2) , h ) if m = 



if m > 



< 



X 



n,(m+l) 



^■n,(m+l) 
L i(n-l),0 



-£(n— l),m 



(31) 



for the 3- and 4-value methods respectively. 



V. PERFORMANCE OF THE ALGORITHMS 



In this section we report on our evaluation of the finite difference algori thm s in 
the previous section, our evaluation criteria being briefly outlined in Section III . We 
first present the results of the simple evaluations checking conservation of energy 
and ergodicity, before presenting more extensive results on phase separation and 
surface tension. 
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A. Energy conservation 



A system consisting of a single simple harmonic oscillator was allowed to evolve 
for a thousand timesteps from an initial configuration with unit position (x = 1). 
Two sizes of timestep were considered: one tenth and one hundredth the period of 
the oscillator. The variation of the energy E = (x 2 + x 2 )/2 relative to the initial 
energy E = | for each of these sets of simulations is shown in Tables | and |l| 
respectively; results are identical for A — \ and A = 1. Where a median and 
standard deviation are given, the energy varies periodically; where an exponent 
is given, the energy grows or decays as aexp(fri). From these tables we can see 
that Gear-3 (q — 1) and Runge-Kutta have an undesirable continual energy loss 
for both sizes of timestep. The other Gear methods show a periodic variation of 
energy for the small timestep, but display a long-term energy drift for the large 
timestep. Both Ferrario methods, Groot- Warren, and Hoogerbrugge-Koelman only 
show periodic variation in energy, with Hoogerbrugge-Koelman showing by far the 
largest magnitude of variation. However, none of the other methods include the 
correct (unit) relative energy within the range of their variation, and so we cannot 
conclude that one periodically varying method is significantly better than any other. 



B. Ergodicity 

Two-dimensional DPD systems of 1600 identical particles were allowed to evolve 
to equilibrium (roughly t = 50), and then the temperature fc^T and pressure P 
were observed for the same length of time again. This was repeated three times for 
each situation from different random initial configurations, and the results averaged. 



The exact parameters used in these simulations are shown in Table III. The units 
of these parameters are those natural to the DPD simulation; although no one has 
yet related them to physical scales of length or time, it is possible to do so if one 
desires to apply the simulation technique to describe a real fluid. These particular 
values were suggested by Hoogerbrugge and Koelman. 

From statistical mechanics, we know that the instantaneous pressure (P) of a 
system in terms of the internal virial (W) and the instantaneous temperature (T) 
are given by 

P = ^Sl^, (32) 

where 

Nk B T=)-Y,m l \± l \ 2 , (33) 



and 



( 34 ) 



W 

2 ^ 

i j>i 

where V is the physical volume of space, = — Xj, and | • | indicates vector 
magnitude. The thermodynamic pressure and temperature are the time averages 
of the instantaneous quantities. Because of the work of Espaiiol, Warren, and 
Coveneyjlacl we can be sure that a temperature exists and is meaningful, at least 
in the continuous-time limit. The pressure was calculated immediately after the final 
force evaluation in each algorithm; the temperature was further evaluated using the 
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final values of velocity. At the end of each simulation, the velocity distribution in the 
x-direction was examined; it was statistically indistinguishable from the Maxwellian 
distribution whenever the temperature converged. 

With the exception of Runge-Kutta, all the methods converge to correct tem- 
perature and pressure as we decrease the timestep. We believe this error is due 
to the unusual timestep-dependent nature of the force, because the Runge-Kutta 
algorithm divides the motion into two steps of two-thirds and then one-third the 
total timestep. Of all the algorithms we consider, only the Runge-Kutta scheme 
subdivides the timestep. The relative proportion of the random component of the 
force changes when we reduce the timestep size, and the algorithm does not properly 
take this into account. 

As we can see in Figs, [l] and (in which is plotted relative to its theoretical 
value), the Groot- Warren (A = ^) and Gear-3 (A = 4, q = 1) methods give the most 
accurate temperature. Figs. || and [| show that the pressure (in DPD units) is most 
accurate for a given timestep with the Gear-3 (A = 1, q = 2) and Gear-4 (A = | 
and A = 1, q = 2) methods, with Gear-3 (A = |, q = 1 and q = 2), Groot- Warren 
(A = |), and Hoogerbrugge-Koelman following close behind. (Groot- Warren (A = 
i) and Hoogerbrugge-Koelman gave nearly identical results for pressure, so their 
symbols are superimposed in Figs. || and [| — look under the curve with the filled 
squares: Gear-3 (A = i, q = 1).) The algorithms are unstable for timesteps larger 
than those shown in Figs. |l| and [| 

Excluding the time taken to calculate the temperature and pressure, these sim- 
ulations take 140 ms per timestep for the Hoogerbrugge-Koelman algorithm on a 
133 MHz DEC Alpha. Because of the large number of particles involved, the speed 
of the algorithms is almost directly proportional to the number of force evaluations 
per timestep: one for Hoogerbrugge-Koelman, the Ferrario methods, and Groot- 
Warren; two for Runge-Kutta; and q for the Gear predictor-correctors. 



C. Phase separation kinetics 



The study of growth kinetics in binary immiscible fluids has received much atten- 
tion lately (see oip.p*pvious papersE3EJ for detailed references, but note also some 
more recent workEJ'Lil) . A central quantity is the characteristic domain size R(t); 
typically, one finds that 

R(t) ~ t 13 . (35) 

Without hydrodynamic interactions, theory and experimental tell us that the scal- 
ing exponent /3 = A. If flow effects are relevant 



( § for R^R h 

(3 = i (in two dimensions) (36) 

1 | for R > R h 



3 



for early-time R <C Rd 
(3 = ^ 1 for late-time Rd <C R <C Rh (in three dimensions) , (37) 
for R > R h 



where Rh = rj 2 /(pr) is the hydrodynamic length and Rd = \fr\D is the diffusive 
length, expressed in terms of the absolute (a.k.a. dynamic) viscosity r\ (dimen- 
sional apaljjsis suggests it is not the kinematic viscosity v = rj/p, as previously 
thoughtc3o), density p, surface tension r, and diffusion coefficient D. 
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Because of the compute-intensive nature of the simulations in this and the follow- 
ing section, we had time enough only , to use a single finite-difference method. Our 
previously performed simulationsE3o used the Hoogerbrugge-Koelman algorithm, 
and of all the methods Groot- Warren (A = h) and Gear-3 (A = h, q = 1) gave 
the most accurate temperature and pressure. We chose to use the Groot- Warren 
algorithm over Gear-3 because of its superior energy conservation. 

Two dimensional systems of 40000 particles were allowed to evolve from a sym- 
metric quench, the initial random configuration differing for each of several simu- 
lations. The model parameters are identical to those used for the smaller simula- 
tions (sec Table III). Observing the growth of those simulations which used the 
Groot- Warren algorithm (A = h, h = 0.1) indicated two regimes of growth fol- 
lowing Eq. ©: (3 = \ (0.478 ± 0.004) crossing over to (3 = § (0.65 ± 0.02) at 
R = 8.7 ± 0.4. Reducing the timestep to h — 0.05, we observe nearly identical re- 
sults: = \ (0.476 ±0.004) crossing to (3 = § (0.63 ±0.03) at R = 7.4 ±0.9. Error 
ranges given are the 68% confidence intervals for the mean, using approximately five 
simulations for each configuration. Phase separation did not occur in symmetric 
quenches with a larger timestep (h — 0.5) for which the algorithm is unstable. This 
is not surprising, as the temperature is undefined and diverges to infinity, so that 
the system does not remain quenched. At this size of timestep, the fluid modeled 
approximates a gas, which we expect to remain in a thoroughly mixed state. 

Asymmetric quenches with a 60:40 majority: minority phase (color) ratio (Groot- 
Warren A = |, h = 0.1) displayed growth with (3 — i (0.476 ± 0.002 crossing to 
0.547 ± 0.003 at R = 7.4 ± 0.7). Asymmetric quenches with a 70:30 ratio (Groot- 
Warren A = |, h = 0.1) also gave (3 — \ (0.460 ± 0.003). Special simulations 
were also set up, in which the conservation of momentum was violated by 10-15% 
for each interaction (Groot- Warren A = h, h = 0.1) by adding a random vector 
to the velocity of each particle every timestep. These simulations demonstrated a 
growth exponent of (3 = 0.407 ± 0.006, ceasing growth entirely at R = 7.57 ± 0.05. 
This suggests that momentum conservation is necessary for the viscous mechanism 
of phase separation. The observed growth exponent is close to the theoretically- 
predicted (3 = |. 

A sample log-log plot of the characteristic domain size (R) against time for a 
single symmetric quench (Groot- Warren A = |, h = 0.1) is shown in Fig. |5|. In this 
figure we can see the initial non-algebraic growth as the system settles down from 
its quench, the (3 = \ early-time diffusive scaling regime, and the (3 = | late-time 
viscous scaling regime. The results presented in this section are qualitatively the 
same as these-pbtained in our earlier simulations using the Hoogerbrugge-Koelman 
algorithm.E2lE-3 



D. Domain surface tension 



Phase separation in binary fluids depends, among other things, on the interfacial 
tension which exists between the two immiscible phases. A further important test 
of our algorithms is thus to check on the existence of a surface tension .by confirming 
the validity of Laplace's law using a series of bubble simulations .E3o 

Our procedure is to set up a circular bubble (radius R) of one color phase within 
the other phase, and allow the bubble to reach equilibrium. We then calculate the 
pressure inside (r < 0.7R) and outside (r > 1.3R) of the bubble, using Eqs. (|32| HJ). 
Repeating these experiments for various size bubbles, we can verify Laplace's law 

T 

■Pin — Pout = (38) 



12 



From our results in Fig. || (Groot- Warren, A = 5, h = 0.1) we can see that we 
have the desired linear behavior, and estimate r, the interfacial surface tension, to 
be 0.33 ± 0.Q5._This is larger than observed in our earlier Hoogerbrugge-Koelman 
simulations ,E3'Eil but because of the greatly decreased noise with the Groot- Warren 
(A = i) method, we are more confident in the accuracy of the current results. 



E. Viscosity 

The absolute (dynamic) viscosity, of this system can be estimated theoretically 
from the continuous-time viscosityEj as 77 = 2.8 ± 0.4, where the error is assumed to 
be of the same order of magnitude as that of the temperature. In order to verify this 
estimate, we performed a series of simulations of steady shear of an homogeneous 
fluid, using Lees-Edwards periodic boundary conditions. Because we intend to 
use this value of viscosity to determine the hydrodynamic length Rh for comparison 



with the spinodal decomposition simulations in Section V C, we continue to use only 
the Groot- Warren algorithm (A = 5). A total of 63 simulations were performed, 
each from a different random initial configuration. Systems of both 1600 and 6400 
particles were studied, six simulations at each of nine different shear rates for the 
former and three simulations at each of three distinct shear rates for the latter. 
As the results from the larger simulations gave a mean viscosity nearly identical 
to that of the smaller ones, we can conclude that domain size effects do not bias 
the smaller, faster simulations. The velocity profile was calculated for each set of 
parameters, and was found to be statistically indistinguishable from linear in every 
case. 

Analyzing these simulations led to a conclusion of r\ = 1.94±0.01. Fig. |?] displays 
the results; the error bars are the 68% confidence intervals for the mean, equally 
weighting each of the simulations for each set of parameters. Others have also found 
discrepancies between theory and simulation, particularly regarding the kinematic 
contribution to viscosityES 

Using the surface tension and viscosity calculated from our simulations, we can 
estimate the hydrodynamic length Rh = rf /(pr) to be 2.9 ± 0.5. During our sim- 
ulations we observed the typical domain size of the crossover in spinodal decompo- 



sition, which we estimate as 8 ± 1 (see Section VC). This length is in agreement 
with Eq. (|36"1), since for length scales much less than Rh we typically see a growth 
exponent of j3 = |, while for length scales much larger than Rh we typically see 
f3 = |. Although it would be more reassuring if Rh and our observed crossover 
length were more comparable in size, we puist remember that the theory leading 
to the proposal of these growth exponentscj makes assumptions about the under- 
lying growth dynamics which break down in this region. This leaves us without an 
accurate theoretical estimate of the characteristic domain size at the crossover. 



F. Radial pair-correlation 

It has been observed that DPD simulations of an ideal gas (i.e. a = 0) using the 
Hoogerbr-ugge-Koelman scheme have unusual structure in the radiaL-pair-correlation 
function.Ej Since continuous-time DPD satisfies detailed balancejlj we expect the 
radial pair-correlation function to have a constant unit value. This suggests that a 
good finite-difference algorithm would not display this unusual structure. 

In order to test each of the finite-difference methods under consideration, systems 
of 40000 identical particles were evolved (a — 0) from initially random configura- 
tions until they reached equilibrium (roughly t = 50), at which point the radial 
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pair-correlation function was calculated. Simulations were performed for each of 
the finite-difference algorithms, using a timestep of h = 0.05 for all methods except 
Ferrario I and II, and Gear-4 (A = 1, q = 1). These results are shown in Fig. pL 
Fig. ^ shows the results of the simulations using the Ferrario methods and Gear-4 
with A = 1 and q = 1, for which we used h — 0.01 since they are unstable at 
the larger timestep. For comparison, the Gear-3 (A = \, q — 1), Groot- Warren 
(A = and Hoogerbrugge-Koelman methods were also tested with the smaller 
timestep (h = 0.01) and so are also shown in Fig. 0. From these figures we can see 
that none of the methods are significantly better than the Hoogerbrugge-Koelman 
scheme, and three are noticeably worse: Gear-3 (A = 1, q = 1), Gear-4 (A = i, 
q = 1), and Groot-Warren (A = 1). Decreasing the timestep gives a marked im- 
provement, so that at h = 0.01 all the methods capture the expected unit radial 
pair-correlation function to a good approximation. 



VI. CONCLUSIONS 



The most significant difference.|-b|etween the simulations reported in this paper 
and those performed previouslyEncj is the improved precision of the measured be- 
havior. We believe this is largely due to changing the discrete-time algorithm from 
Hoogerbrugge-Koelman to Groot-Warren. Coupled with its superior performance 
in the less-rigorous tests, we recommend the Groot-Warren method with A = ^ and 
timestep h < 0.1 as most suitable for DPD simulations. The Gear-3 algorithm with 
A = I and q = 1 is the obvious second choice, as it gave temperature and pressure 
to nearly the same accuracy as the Groot-Warren method (A — i), and is stable to 
a timestep nearly as large. Our only concern in recommending this modified Gear 
algorithm is its poor energy conservation. 

We suspect the main limitation to improving the existing algorithms is their in- 
ability to properly address the stochastic component of the forces; there is some 
literature on Brownian dynamics simulationsEj that may be useful in finding an 
even better algorithm. Compared with the solution of ordinary differential equa- 
tions, the numerical solution of stochastic differential equations-.is a recent area of 
mathematics, and one in which much more research is needed.Ej 
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Method 






Median 


Std. dev. 


Exponent 


Ferrario I (Ec 
Ferrario II (E 
Gear-3 (Eq. ( 
Gear-3 (Eq. ( 
Gear-4 (Eq. ( 
Gear-4 (Eq. ( 
Groot-Warrer 
Hoogerbrugge 
Runge-Kutta 
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50 

'£ 

31 

~ 

-b 

(I 


M) 

),9 = 1) 
),9 = 2) 
),9 = 1) 

Eq. ©) 

Coelman (Eq. (gj)) 
]q. ©) 


1.142 
1.054 

1.054 
1.107 


0.085 
0.039 

0.039 
0.25 


-0.024 
0.0016 
-0.0044 
0.00017 

-0.0044 


TABLE I. Simple harmonic oscillator energy relative to initial energy for an homo- 
geneous fluid, with a timestep of one tenth the period. A given median and standard 
deviation indicate periodic variation; an exponent describes exponential growth or decay. 
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Median 


Std. dev. 


Exponent 


Ferrario I (Ec 
Ferrario II (E 
Gear-3 (Eq. ( 
Gear-3 (Eq. ( 
Gear-4 (Eq. ( 
Gear-4 (Eq. ( 
Groot-Warrer 
Hoogerbrugge 
Runge-Kutta 


1- 
q. 

5C 

n 

31 

'( 
-b 

(1 


@) 
If) 

),9 = 1) 
), 9 = 2) 
), 9 = 1) 

Eq. ©) 

Coelman (Eq. (gj)) 
]q. ©) 


1.00050 
1.00049 

0.99951 
1.00032 
1.00033 
1.00049 
1.0020 


0.00042 
0.00035 

0.00012 
0.000012 
0.000013 
0.00035 
0.022 


-2.6 x 10 -6 
-4.4 x 10" 7 



TABLE II. Simple harmonic oscillator energy relative to initial energy for an homoge- 
neous fluid, with a timestep of one hundredth the period. A given median and standard 
deviation indicate periodic variation; an exponent describes exponential decay. 





Model 




Parameter 


Value 


QO 


7.063 


ai 


7.487 


7 


5.650 


a 


1.290 


in, 


1 


r c 


1.3 


P 


4 



TABLE III. DPD model parameters (see Section |[l|). 
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FIG. 1. Relative equilibrium temperature (fcsT) vs. timestep (h) for an homogeneous 

fluid. Data are shown for the various methods: 1 Ferrario I, — x — Ferrario II, 

------- Gear-3 (A = 1, q = 1), Gear-3 (A = i, q = 1), ---o--- Gear-3 

(A = 1, q = 2), — •— Gear-3 (A = i, q = 2), — A — Gear-4 (A = 1, q = 1), 
---A--- Gear-4 (A = ±, q = 1), • • • V • • • Gear-4 (A = 1, q = 2), Gear-4 (A = ±, 

q — 2), — — Groot-Warren (A = 1), Groot- Warren (A = |), — * — Hooger- 

brugge-Koelman, ■ ■ -k ■ ■ Runge-Kutta. The solid line drawn at unity is for reference. 



FIG. 2. Relative equilibrium temperature (fcsT) vs. timestep (h) for an homogeneous 
fluid (detail). Data are shown for the various methods: — X-- Ferrario II, — □ — Gear-3 
(A = 1, q = 1), Gear-3 (A = ±, q = 1), — A — Gear-4 (A = 1, q = 1), 

---A--- Gear-4 (A = ±, q = 1), —<)— Groot-Warren (A = 1), — ♦ — Groot-Warren 
(A = |). The solid line drawn at unity is for reference. 

FIG. 3. Equilibrium pressure (P) vs. timestep (h) for an homogeneous fluid. Data are 

shown for the various methods: 1 Ferrario I, — x — Ferrario II, — □ — Gear-3 

(A = 1, q = 1), Gear-3 (A = \, q = 1), ---o--- Gear-3 (A = 1, q = 2), 

— •— Gear-3 (A = ±, q = 2), — A — Gear-4 (A = 1, q = 1), ---A--- Gear-4 

(A = i, q = 1), ...v--- Gear-4 (A = 1, q = 2), T Gear-4 (A = ±, 

q = 2), — — Groot-Warren (A = 1), Groot-Warren (A = i), — * — Hooger- 

brugge-Koelman, + Runge-Kutta. 



FIG. 4. Equilibrium pressure (P) vs. timestep (h) for an homogeneous fluid (detail). 

Data are shown for the various methods: Gear-3 (A = |, q = 1), o Gear-3 

(A = 1, q = 2), — •— Gear-3 (A = ±, q = 2), ..-v--- Gear-4 (A = 1, q = 2), 

T Gear-4 (A = |, q = 2), — ♦-- Groot-Warren (A = |), ---*--- Hooger- 

brugge-Koelman, Runge-Kutta. 



FIG. 5. Growth of characteristic domain size for a single symmetric quench 
(Groot-Warren A = \, h = 0.1). Solid lines are of slope | and |. 



FIG. 6. Pressure difference as a function of bubble size (Groot-Warren, A = |, h = 0.1). 



FIG. 7. Absolute viscosity (r/) vs. shear rate for an homogeneous fluid. Symbols indicate 
simulation size: + N = 1600, and x N = 6400. 



FIG. 8. Radial pair correlation function g(r) vs. distance r for a DPD ideal gas (a = 0), 
with timestep h = 0.05. Data are shown for the various methods: — □ — Gear-3 (A = 1, 
q = 1), Gear-3 (A = ±, q = 1), -.-o--- Gear-3 (A = 1, q = 2), — •— Gear-3 

(A = i, q = 2), ---A--- Gear-4 (A = ±, q = 1), • • • V • • • Gear-4 (A = 1, q = 2), 

T Gear-4 (A = \, q = 2), —<>— Groot-Warren (A = 1), — ♦ — Groot-Warren 

(A = i), — * — Hoogerbrugge-Koelman, Runge-Kutta. The solid line drawn at 

unity is for reference. 
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FIG. 9. Radial pair correlation function g(r) vs. distance r for a DPD ideal gas (a = 0), 

with timestep h = 0.01. Data are shown for the various methods: 1 Ferrario I, 

— x — Ferrario II, Gear-3 (A = ^, q = 1), — A-- Gear-4 (A = 1, q = 1), 

Groot-Warren (A = i), — * — Hoogerbrugge-Koelman. The solid line drawn at 
unity is for reference. 
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FIG. 1. Keir E. Novik, Journal of Chemical Physics 
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FIG. 2. Keir E. Novik, Journal of Chemical Physics 
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FIG. 3. Keir E. Novik, Journal of Chemical Physics 
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FIG. 4. Keir E. Novik, Journal of Chemical Physics 
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FIG. 5. Keir E. Novik, Journal of Chemical Physics 
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FIG. 6. Keir E. Novik, Journal of Chemical Physics 
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FIG. 7. Keir E. Novik, Journal of Chemical Physics 
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FIG. 8. Keir E. Novik, Journal of Chemical Physics 
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FIG. 9. Keir E. Novik, Journal of Chemical Physics 
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